function value = f(p)
x = p(:,1); y = p(:,2);
value = (x.^2 + y.^2).*(1./((x - 1/2).^2 + (y - 1/2).^2 + 1/100)...
        - 1./((x + 1/2).^2 + (y + 1/2).^2 + 1/100))...
        - 10.*sin(y).*((2.*y - 1)./((x - 1/2).^2 + (y - 1/2).^2 + 1/100).^2 ...
        - (2.*y + 1)./((x + 1/2).^2 + (y + 1/2).^2 + 1/100).^2)...
        - 10.*cos(y).*((2.*(2.*x - 1).^2)./((x - 1/2).^2 + (y - 1/2).^2 ...
        + 1/100).^3 - (2.*(2.*x + 1).^2)./((x + 1/2).^2 + (y + 1/2).^2 ...
        + 1/100).^3 - 2./((x - 1/2).^2 + (y - 1/2).^2 + 1/100).^2 ...
        + 2./((x + 1/2).^2 + (y + 1/2).^2 + 1/100).^2) ...
        - 10.*cos(y).*((2.*(2.*y - 1).^2)./((x - 1/2).^2 ...
        + (y - 1/2).^2 + 1/100).^3 - (2.*(2.*y + 1).^2)...
        ./((x + 1/2).^2 + (y + 1/2).^2 + 1/100).^3 - 2./((x - 1/2).^2 ...
        + (y - 1/2).^2 + 1/100).^2 + 2./((x + 1/2).^2 + (y + 1/2).^2 ...
        + 1/100).^2);
end 